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I.  INTRODUCTION 


Problems  involving  solidification  or  melting  of  materials  are  of  con¬ 
siderable  importance  in  many  technical  fields  such  as  casting  processes, 
food  freezing,  crystal  growing,  cryosurgery,  thermal  energy  storage,  and 
satellite  temperature  control  applications,  to  mention  a  few.  Heat  conduction 
problems  involving  freezing  or  melting  are  complicated  due  to  the  coupling 
of  the  temperature  field  with  the  rate  of  propagation  of  the  phase  boundary 
between  the  solid  and  liquid  phases.  Only  a  few  exact  analytical  solutions 
have  been  found,  e.g.,  Neumann's  solution  to  the  one-dimensional  Cartesian 
case  reported  in  Carslaw  and  Jaeger.^  Most  available  solutions  are  obtained 
by  analytical  approximations  and  numerical  methods.  The  one -dimensional 
problem  for  simple  shapes  such  as  plates,  cylinders,  and  spheres  has  been 

treated  by  various  analytical  techniques,  e.g.,  the  heat  balance  integral 

2  3 

method  by  Goodman,  the  variational  method  by  Biot,  the  method  of  moving 

4 

heat  sources  by  Rosenthal  and  the  method  of  polynomial  approximation  by 
5.6  7 

Megerlin  and  Lin.  London  and  Seban  obtained  exact  closed-form  solutions 
by  assuming  that  the  heat  capacity  of  the  solidified  or  melted  substance  is 
negligible  relative  to  the  latent  heat  of  fusion.  This  assumption  greatly 


Carslaw,  H.S.,  and  Jaeger,  J.  C.,  Conduction  of  Heat  in  Solids,  2nd  ed., 
Oxford  University  Press,  London  and  New  York,  1959. 

2 

Goodman,  T.  R.,  "The  Heat-Balance  Integral  and  Its  Application  to  Problems 
Involving  a  Change  of  Phase,  "  Transactions  of  ASME,  Vol.  80,  1958, 
pp.  335-342. 

3 

Biot,  M.  A.,  and  Daughaday,  H.,  "Variational  Analysis  of  Ablation,  " 

Journal  of  Aerospace  Sciences,  Vol.  29,  No.  2,  1962,  pp.  227-228. 

4 

Rosenthal,  D.,  "The  Theory  of  Moving  Sources  of  Heat  and  Its  Application 
to  Metal  Treatments,  ”  Transactions  of  ASME,  Vol.  68,  1946,  pp.  849-866. 

5 

Megerlin,  F.,  "Geometrisch  eindimensionale  Warmeleitung  beim  Schmelzen 
und  Erstarren,  "  Forsch.  Ing.-Wes.,  Vol.  34,  1968,  pp.  40-46. 

Lin,  S.,  "An  Analytical  Method  for  Solving  Geometric  One-dimensional 
Freezing  or  Melting  Problems,  "  ASME  Paper  No.  73-WA/HT-33. 

7 

London,  A.  L.,  and  Seban,  R.A.,  "Rate  of  Ice  Formation,  "  Transactions  of 
the  ASME,  Vol.  6,  1943,  pp.  771-778. 
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simplifies  the  matheinatics  but  limits  the  applicability  of  the  solution  to  a 
certain  class  of  problems  as  discussed  below.  Numerical  solutions  to  the 
one -dimensional  phase-change  problems  for  simple  shapes  such  as  slabs, 
cylinders,  and  spheres  are  too  numerous  to  mention  here.  References  8  and 
9  are  good  starting  points. 

Multidimensional  phase-change  problems  are  very  complex  to  solve 

analytically  even  for  highly  idealized  situations  such  as  those  treated  in 

1 2 

Refs.  10  and  11;  the  common  approach  is  to  use  numerical  methods.  How¬ 
ever,  several  seemingly  multidimensional  problems  can  be  closely  approxi¬ 
mated  by  one-dimensional  techniques  if  a  proper  coordinate  system  is  used  so 

that  the  physical  boundaries  of  the  system  can  be  made  to  coincide  with  the 

13  14 

coordinate  surfaces.  ’  The  purpose  of  this  report  is  to  develop  a  unified 
approach  in  the  use  of  orthogonal  curvilinear  coordinates  to  conduct  prob¬ 
lems  involving  phase-change  and  to  present  examples  in  several  of  the  most 
common  coordinate  systems.  The  closed-form  expressions  obtained  give 
the  solidified  (or  melted)  fraction,  interface  position,  boundary  temperatures, 
and  heat  fluxes  as  a  function  of  time  for  elliptic  cylinder  and  spheroidal  phase - 
change  containers  of  various  eccentricities. 


g 

Ehrlich,  L.  W.,  "A  Numerical  Method  of  Solving  a  Heat  Flow  Problem  with 
Moving  Boundary,  11  Journal  of  ACM,  Vol.  5,  1958,  pp.  161-176. 

9 

Murray,  W.  D.,  and  Landis,  F.,  "Numerical  and  Machine  Solutions  of 
Transient  Heat -Conduction  Problems  Involving  Melting  or  Freezing,  " 
Transaction  of  ASME,  Vol.  81,  1959,  pp.  106-112. 

Poots,  G.,  "An  Approximate  Treatment  of  Heat  Conduction  Problem  In¬ 
volving  a  Two-Dimensional  Solidification  Front,  "  International  Journal  of 
Heat  and  Mass  Transfer,  Vol.  5,  1962,  pp.  339-348. 

Rathjen,  K.A.,  and  Jiji,  L.  M.,  "Heat  Conduction  with  Melting  or  Freezing 
in  a  Corner,  "  Journal  of  Heat  Transfer,  Trans.  ASME,  Series  C.  Vol.  93 
1971,  pp.  101-109. 

2 

Shamsundar,  N.,  and  Sparrow,  E.  M.,  "Analysis  of  Multidimensional 
Conduction  Phase  Change  Via  the  Enthalpy  Model,  "  Journal  of  Heat  Transfer 
Trans.  ASME.  Series  C,  Vol.  97,  1975,  pp.  333-3407"  * 

3 

Yovanovich,  M.  M.,  Advanced  Heat  Conduction,  Hemisphere  Publishing 
Corporation,  Washington,  D.C.,  1978. 

4 

Moon,  P.,  and  Spencer,  D.  E.,  Field  Theory  for  Engineers.  D.  Van  Nostrand 
Company,  Inc.,  Princeton,  New  Jersey,  1961. 
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II.  ANALYSIS 


A.  CARTESIAN  COORDINATES 

Consider  first  the  relatively  straightforward  case  of  idealized  freezing 
between  two  infinite  parallel  plates  as  shown  in  Fig.  la.  At  time  t  =  0  the 
content  between  parallel  plates  is  assumed  to  be  a  pure  liquid  at  its  freezing 
or  melting  temperature.  As  the  outer  walls  are  cooled  by  contact  with 
another  medium  at  T  ,  the  solidification  boundary  moves  inward  toward 
smaller  values  of  x.  The  effects  of  supercooling,  density  changes  due  to 
change  of  phase,  convection  between  the  two  phases,  and  temperature  drop 
through  the  container  wall  are  neglected.  Let  Q  be  the  heat  rate  removed 
from  the  surface  at  x  through  a  heat  transfer  coefficient  hQ.  If  the  sensible 
thermal  energy  stored  in  the  solidified  layer  is  small  compared  to  the  latent  - 
heat-of-freezing,  heat  flowing  from  the  solid-liquid  interface  at  x  is  also 
equal  to  Q  and  is  given  by 


o  o 


This  heat  flow  provides  extraction  of  the  latent  heat  of  fusion  necessary  for 
freezing  at  the  surface  x 

Q  =  L^  (xo  -  x)A p  (2) 

where  (d/dt)(x  -  x)A p  is  the  mass  rate  of  solidification  at  the  growing  surface 
kg/s,  and  L  is  the  latent  heat  of  fusion  cal/gm.  When  Eqs.  (1)  and  (2)  are 
equated  to  eliminate  Q 
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where  A  =  Aq  for  Cartesian  coordinates.  Equation  (3)  can  be  readily  non - 
dimensionalized  by  employing  the  dimensionless  distance  coordinate 

>Jj 

x  =  x/x  £1.0  and  the  Biot  number  Bi  =  h  x  /k  as  follows: 
o  o  o 


k(T,  -  T  ) 
f  r _ o 


dt 


(4) 


For  the  case  of  inward  solidification,  Eq.  (4)  can  be  integrated  from  x  =  1 
at  t  =  0  to  x  =x  at  t;  the  result  is 


k(Tf  -  T  )t 
fr  o 


=  SteFo 


(5) 


where  the  dimensionless  time  parameter  is  the  product  of  the  Stefan  number 

Ste  -  Tq)/L  times  the  Fourier  number  Fo  =  (k/pCL  )t.  The  Stefan 

number  is  the  ratio  of  the  specific  heat  effect  to  Ktent  heat.  For  heat  storage 

applications  it  is  sometimes  more  convenient  to  express  Eq.  (5)  in  terms  of 

the  frozen  fraction  F,  which  is  defined  by  F  =  (x  -  x)/x  =  1  -  x 

7  o  o 


Ste  Fo 


(6) 


The  wall  temperature  history  can  be  readily  obtained  by  equating  the  heat 
flow  through  the  outer  surface  hQAQ(Tw  -  Tq)  to  Eq.  (1)  and  rearranging 


1 _ 

1  +  Bi  (1  -  x*) 


1 

1  +  BiF 


/T 


2 

+  2Bi  SteFo 


(7) 
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The  plotted  results  from  Eqs.  (6)  and  (7)  will  be  presented  later  when 
freezing  in  elliptic  cylinders  is  compared  with  oblate  spheroids  of  high  eccen¬ 
tricity.  Note  that  the  boundary  condition  temperature  Tq  in  Eq.  (2)  could 
be  a  function  of  time  without  complicating  the  integration  significantly. 
Closed-form  solutions  are  still  obtainable,  as  long  as  the  assumed  function 
can  be  integrated. 

Equations  (6)  and  (7)  are  also  applicable  for  the  case  of  melting,  as 

shown  in  Fig.  lc,  the  only  difference  being  that  in  this  case  F  is  the  melted 

fraction  and  T  >  T  >  T,  .  The  case  of  freezing  in  an  outward  direction, 
o  w  fr 

Fig.  lb,  can  also  be  similarly  derived  by  equating  the  heat  extraction  rate 
to  the  mass  rate  of  freezing 


T 


fr 


-  T 

o 


1 

h  A 


o 


_  .  t  dx 

pALar 


(8) 


Using  the  same  nondimensional  terms  as  before  and  integrating  from  x  =  1 
at  t  =  0  to  x  =  x  at  t  results  in 

SteFo  =  (  -j=p- 

\Bl 


where  x  =x/xq  ^  1.0.  This  equation  is  valid  for  outward  freezing  or  melting 
(see  Figs,  lb  and  Id). 


B.  GENERALIZED  COORDINATES 

Before  proceeding  with  formulation  of  the  problem  in  general  curvilinea 
coordinates,  a  general  expression  for  the  thermal  resistance,  comparable  to 
the  x/kA  for  Cartesian  coordinates,  is  required.  Special  coordinate  systems 
are  used  because  it  is  sometimes  possible  to  make  the  isothermal  surfaces 
of  a  given  heat  conduction  problem  coincident  with  coordinate  surfaces,  thus 
making  the  temperature  field  one -dimensiona  1  in  that  coordinate  system,  e.g. 
the  field  between  two  concentric  cylinders  in  circular  cylinder  coordinates. 
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Consider  the  orthogonal  curvilinear  coordinates  shown  in  Fig.  2.  The 
general  line  element  ds  is  the  diagonal  of  the  infinitesimal  parallelepiped 
with  faces  that  coincide  with  the  planes  Uj,  u2,  or  =  const  and  is  given  by 

ds2  =  gjdu2  + g2du2  + g3du2  (10) 

which  can  be  considered  as  the  definition  of  the  metric  coefficients  g2.g2.g3 
that  may  be  functions  of  u^,  u2>  and  u^.  When  the  Cartesian  coordinates 
x,  y,  z  can  be  expressed  in  terms  of  the  new  coordinates,  u^,  u2,  u^  by  the 
equations  x  =  x(u2>  u2,  u^),  y  =  yfu^,  u2,  u^),  and  z  =  z (u^ ,  u2>  u^ ),  the  metric 
coefficients  can  be  readily  generated  by  means  of  the  following  formula:^2'  ^ 


gi  =  dr)  +(t£r)  +  (It) 


i  =  1,2,3 


The  infinitesimal  volume  is  given  by 

dV  =  \fg  du^du2du3  (12) 

where  sfg  =  y/g^g^g^" .  The  elemental  surface  area  orthogonal  to  the 
u^-direction  is 

dAj  =  y  g2g3  du2du3  (13) 

Simila-r  expressions  can  be  written  for  the  areas  in  the  u2~  and  ^-directions. 
The  heat  flow/unit  time  through  this  surface  into  the  volume  element  is 


dQ.  =  -  kdA.  ^ 
1  Ids, 


,  >Jg  dT  ,  , 

-  k  — -T —  du^du, 
gi  du.  23 


The  net  rate  of  heat  conduction  out  of  the  volume  element  is 


JL.  /k^L  dT_\ 

dul  \  gl  dul/ 


dujdu2du3 


(15) 


z 


Fig.  2.  Orthogonal  Curvilinear  Coordinates 


Dividing  by  volume  element  dV  =  \fg  du^du^du^  and  equating  to  zero  results  in 
Laplace's  equation  in  the  u^  direction 


J_  _£L  L&  £l)  -  o 

/g  d“l  y  «1  *»l/ 


Integrating  once,  the  result  is 

k  ^  =  C.  (] 

gi  duj  1 

At  this  point  the  thermal  conductivity  is  assumed  to  be  constant  and 
the  boundary  conditions  are  taken  as 


T  =  T, 


T  =  T. 


ul  =  a 


Uj  =  b 


Integrating  Eq.  (17)  between  these  limits  and  solving  for  results  in 


C1  = 


k(T2  -  T x ) 


Substituting  in  Eq.  (17)  and  solving  for  the  temperature  gradient  results  in 


T2-Tl  *1 


/S*. 


The  heat  flow/unit  time  along  the  u^ -direction  can  be  obtained  by  substituting 
in  Eq.  (14)  and  integrating  between  appropriate  limits  along  u^  and  u^ 


Ql  =  k(Tl  -  T2) 


j  J  du2du3 

"2  /  5s 
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The  thermal  resistance  R  is  defined  as  the  temperature  drop  over  the 
total  heat  flow  rate;  thus 


(22) 


v  ,  13, 15  , 

Yovanovich  derived  this  equation  for  the  more  general  case  of  thermal 

conductivity  as  a  linear  function  of  temperature  and  showed  that  Eq.  (22)  is 
valid  if  k  represents  the  arithmetic  average  conductivity.  He  also  specialized 


Eq.  (22)  for  the  most  frequently  used  coordinate  systems 
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Consider  a  phase-change  material  (PCM)  confined  between  two  iso¬ 
thermal  surfaces  of  a  general  orthogonal  coordinate  system  at  different 
temperature  levels  with  phase  change  taking  place  inwardly  from  a  larger 
outer  surface  to  a  smaller  inner  surface.  A  heat  balance  similar  to  Eq.  (3) 
for  the  Cartesian  coordinate  case  can  be  written  as  follows: 


where  all  assumptions  made  in  deriving  Eq.  (3)  apply  for  this  equation  as  well. 
When  phase  change  proceeds  in  an  outward  direction,  from  the  smaller  area 
toward  a  larger  area,  the  right-hand  side  of  the  equation  is  changed  to  pLdV/dt. 
Equation  (23)  will  be  specialized  and  solved  for  elliptic  cylinder,  oblate 
spheroidal,  prolate  spheroidal,  and  bicylindrical  coordinates. 

1  5Yovanovich,  M.M.,  "A  General  Expression  for  Predicting  Conduction 
Shape  Factors,  "  AIAA  Paper  No.  73-121. 
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c.  ELLIPTIC  CYLINDER  COORDINATES  (?7,  fl,  z) 

Phase  change  conduction  problems  in  or  around  elliptic  cylinders  can  be 
enormously  simplified  by  using  elliptic  cylinder  coordinates  as  shown  in 
Fig.  2.  The  coordinates  are  obtained  by  taking  an  orthogonal  family  of 
confocal  ellipses  and  hyperbolas  in  a  plane  and  translating  them  in  the 
z-direction.  The  coordinate  surfaces  of  interest  here  are  elliptic  cylinders 
(T]  =  const).  The  relations  between  the  elliptic  and  Cartesian  coordinates  are 

x  =  a  cosh/7  cos 

y  =  a  sinh/7  sin *P 

z  =  z  (24) 

The  metric  coefficients  can  be  easily  derived  by  the  use  of  Eq.  (11) 

2  2  2 

8n  =  glf>  ~  a  COsh  11  '  COS  gZ  =  1  (25) 

From  Eq.  (24)  it  can  be  seen  that  b  =  a  cosh/7,  c  =  a  sinh/7  and  TJ  =  tanh_1c/b. 
The  thermal  resistance  between  two  elliptic  cylinders  is  readily  obtained  by 
Eq.  (22)  or  taken  directly  from  Refs.  13  or  15. 

% 

R  “  27r£k  (26> 

The  heat  balance  equation  is 


=  PLdt  (Vo  -  V) 


The  area  Aq  can  be  computed  from  Eq.  (13)  as  follows: 

/  27T  l  7T/2 

Ao  =  f  f  /^^dz  =  4  J  J  a  'xosh^-cos^  Ip  dl/'d  5 

0  0  00 


* 


Lett-ng  ^  -  90  -  0,  it  can  be  shown  that 

Ao  =  4a^  cosh*70  E(»70)  (29) 

where  E(77q)  is  the  complete  elliptic  integral  of  the  first  kind  tabulated  in  most 
mathematical  tables  and  shown  in  Fig.  3  for  ease  of  reference.  The  expres¬ 
sion  for  the  volume  V  can  be  calculated  from  Eq.  (12)  by  triple  integration. 

/  7T/2  rj 

v  =  4&Z  f  f  f  (cosh277  -  cos2^)  drjd'l'dz  =  fa2/  Sinh2  T]  (30) 
o  o  o 

Differentiating  with  respect  to  time  and  substituting  the  heat  balance  equation 
becomes 


where  the  semiminor  axis  c  has  been  used  to  define  the  Biot  number  Bi  =  h  c/k. 
The  term  on  the  left-hand  side  of  the  equation  is  the  product  of  the  Stefan  and 
f  ourier  numbers.  Performing  the  integrations  shown  finally  results  in 


* 
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-  1 


k(T,  -  T  )t  n\  (~) 

steFo  =  - r—^-  =  8B.V^h'  ,si„h2,o  -  •inhZII) 

#o  o 


pLc 


(I)  - 


(ZTj^sinhZTj  -  2 77  sinh277  -  cosh2 77q 


+  cosh277) 

Since  b  =  a  cosh^  and  c  =  a  sinh77,  it  can  be  easily  shown  that 

_  _  1  4  b/c  +  1 
V  "  2  b/c  -  1 


(33) 


(34) 


a  = 


(35) 


The  wall  temperature  history  can  be  obtained  by  equating  the  heat  flow  through 

the  surface  h  A  (T  -  T  )  to  the  left-hand  side  of  Eq.  (27)  and  rearranging; 
o  o  w  o 

the  result  is 


1+2|i(n 

7T  o 


-  rj)coshr)oE(J7o) 


(36) 


This  ratio  can  also  be  interpreted  as  the  normalized  heat  flow  rate  that  can  be 
extracted  from  the  PCM  during  freezing  or  the  heat  flow  rate  that  must  be 
introduced  to  achieve  a  given  percent  of  fusion.  Figure  4  shows  the  variation 
of  the  dimensionless  wall  temperature  and  frozen  fraction  F  as  functions  of 
the  dimensionless  time  parameter  SteFo.  The  frozen  fraction  is  related  to 
the  77-coordinate  by 


sinh2?7 

sinh277 

o 


(37) 
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Fig.  4.  Frozen  Fraction  and  Dimensionless  Wall  Temperature  as 
Function  of  Time  for  Elliptic  Coordinates 


As  c/b  approaches  unity,  the  results  approach  the  limiting  case  of  a  circular 

7 

cylinder  derived  by  London  and  Seban.  In  terms  of  the  symbols  used  in 
this  report,  the  circular  cylinder  results  are 


SteFo 


where  r  =  r/r  =1.0  and  the  frozen  fraction  F  =  1 
o 

wall  temperature  is  given  by 


(38) 

The  dimensionless 


1 

1  -  Bi inr 


(39) 


The  results  for  a  slab,  Eqs.  (6)  and  (7),  are  also  shown.  The  case  for 
c/b  =  0.01,  not  shown,  is  closely  equal  to  the  c/b  =  0.  1  case;  thus,  the 
slab  results  do  not  provide  a  good  approximation  for  elliptic  cylinders  of 
high  eccentricities.  Note  that  for  applications  such  as  energy  storage  or 
close  temperature  control,  where  a  low  temperature  difference  is  required 
throughout  the  PCM,  elliptic  cylinder  containers  provide  a  better  choice  than 
circular  cylinders  as  seen  by  comparing  the  c/b  =  0.  9  and  cylinder  results 
in  Fig.  4.  The  effect  of  Biot  number  on  elliptic  cylinders  of  moderate 
eccentricities  is  summarized  in  Fig.  5.  For  applications  where  low  tempera¬ 
ture  differences  are  required,  Biot  number  should  be  kept  as  low  as  possible 
within  the  constraint  of  being  able  to  achieve  complete  phase  change  in  the 
allocated  time  span.  The  effect  of  eccentricity  c/b  on  the  time  required  for 
complete  solidification  and  the  corresponding  dimensionless  wall  tempera¬ 
ture  is  shown  in  Figs.  6  and  7.  Low  Biot  numbers  and  c/b  ratios  result  in 
small  temperature  differences  during  the  phase-change  process;  however,  the 
corresponding  times  for  complete  phase  change  are  relatively  larger.  When 
the  Biot  number  is  sufficiently  high,  a  relatively  large  temperature  difference 
develops  across  the  PCM  for  all  c/b  values. 
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Frozen  Fraction  and  Dimensionless  Wall  Temperature  as 
Function  of  Time  with  Biot  Number  as  a  Parameter 


ELLIPTIC  CYLINDER 


Fig.  7.  Dimensionless  Wall  Temperature  as  a  Function  of 
c/b  with  Biot  Number  as  a  Parameter 
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For  phase-change  problems  in  an  outward  direction  where  T)  >  rj  it 
can  oe  shown  that  the  dimensionless  time  to  reach  a  given  position  rj  is 
obtained  by  placing  a  negative  sign  on  the  first  term  only  of  Eq.  (33).  The 
wall  temperature  of  the  inner  boundary  can  be  calculated  from  Eq.  (36)  if 
(r}Q  -  /?)  is  changed  to  (>7  -  77q).  The  resulting  expression  can  be  used  to 
estimate  the  rate  at  which  a  phase-change  front  would  move  into  the  half¬ 
space  above  the  x-axis  from  a  strip  of  length  2a  located  on  the  x-axis,  which 
is  assumed  to  be  insulated  everywhere  else. 

D.  OBLATE  SPHEROIDAL  COORDINATES  (/?,#>) 

The  oblate  spheroidal  coordinate  system  is  generated  by  taking  an 

orthogonal  family  of  confocal  ellipses  and  hyperbolas  and  rotating  it  about 

the  minor  axis  of  the  ellipses  as  shown  in  Fig.  2.  The  T]  -  const  solids  are 

called  oblate  spheroids.  When  =  0,  the  spheroid  degenerates  into  a  flat 

disk  of  radius  a.  The  relations  between  oblate  spheroidal  coordinates  and 

13 

Cartesian  coordinates  are 

x  =  a  cosh?;  sindcos^ 
y  =  a  cosh T)  sin 6  sin^ 

z  =  a  sinhrj  cos#  (40) 


The  metric  coefficients  can  be  easily  derived  by  using  Eq.  (11) 


g 

g 


n 


2 

a 


=  a^(cosh ^77  -  sin^#) 
cosh^t)  sin^# 


(41) 


From  Eq.  (40)  the  semimajor  axis  is  b  =  a  coshr;  and  the  semiminor  axis  is 
c  =  a  sinht)  from  which  r]  =  tanh  ^c/b. 


-24- 


The  thermal  resistance  between  two  oblate  spheroids  can  be  obtained  by 
performing  the  integrations  indicated  in  Eq.  (22)  or  by  using  the  results  of 
Ref.  15  directly 


tan  *(sinh?7o)  -  tan  \sinh77) 
47Tka 


where  7?  <  T)  for  the  inward  solidification  case.  The  heat  balance  equation 
can  now  be  written  as  before 


where  the  definite  integral  Iq (77)  is  plotted  in  Fig.  3.  Volume  V  can  be  found 
by  integrating  Eq.  (12) 

3 

4-rra  7 

V  =  — ^ —  (sinh77  +  sinh  77)  (45) 


Differentiating  with  respect  to  time  and  substituting  in  Eq.  (43)  above,  to¬ 
gether  with  Eq.  (44),  gives 


k(T  -  T  )t 
ir  o 

r~2 — 

pLc 


/l  >  v(ff- 

(Bi/ 


cosh  T)  I  (rj  ) 
'o  o  'o 


(!)  - 


tan  ^(sinh77Q) 


/  (co 


sh3rj 


r  2  1  770 

cosh/7)d/7  -  M  k  )  -1  •  J  tan  *(sinh77) 

J^cosh3^  -  j  cosh  77  J  drj 


(46) 


The  integrations  shown  on  the  right-hand  side  of  the  equation  can  be  carried 
out  in  closed  form  with  the  following  final  result: 


SteFo  = 


(sr) 


!) 


-  l 


cosh  T]olo(T]o) 


(!)  - 


tan  ^(sinh?7o) 


(!)  - 


(47) 


whe  re 


Ij  =  ~Y2  (s mh3 J7n  +  sinhr/^  -  sinh3t)  -  sinh77) 


=  Tz  [(si 


-l, 


inh3  77Q  +  sinh77Q)tan  (sinh^) 


(sinh3r)  +  sinht))  •  tan 


-  -jlj-  (cosh277o  -  cosh2»7) 


1  (sinhrj)J 


As  in  the  elliptic  cylinder  case,  the  Biot  number  is  defined  in  terms  of  the 
semiminor  axis  c  which  is  related  to  a  by  a  =  'J b2  -  c2. 
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The  walL  temperature  history  can  be  obtained  by  equating  the  heat 

flow  through  the  outer  surface  h  A  (T  -  T  )  to  the  left-hand  side  of 
6  o  o  w  o 

Eq.  (43)  and  rearranging;  the  result  is 


1  + 


2  r  - 1 

Bicosh  tan  (sinhl^) 


tan  *(sinhT7) 


-1 


(48) 


Figure  8  shows  Eqs.  (47)  and  (48)  plotted  for  Bi  =  1.0  for  oblate  spheroids  of 
c/b  ratios  of  0.  9,  0.5,  and  0.  1.  The  frozen  fraction  F  is  related  to  the 
^-coordinate  by 


F 


1.0 


2 

sinh T)  cosh  T] 
sinh^cosh^^ 


(49) 


As  c/b  approaches  unity,  the  results  approach  the  limiting  case  of  a  sphere 
derived  in  Ref.  7.  In  terms  of  the  symbols  used  herein,  the  sphere  results 
are 


SteFo 


'  (bV  ‘)(‘  -  )  +  2  • 

...3 


where  r  =  r/r  S  1.0  and  the  frozen  fraction  F  =  1  -  r 
o 

less  wall  temperature  is  given  by 


(50) 


The  dimension¬ 


al) 


The  case  for  c/b  =  0.  01,  not  shown,  is  almost  equal  to  the  c/b  =  0.  1  case; 
thus  the  slab  results  do  not  provide  a  good  approximation  for  disclike  oblate 
spheroids  of  high  eccentricities.  For  applications  where  a  low  temperature 


Frozen  Fraction  and  Dimensionless  Wall  Temperature  as 
Function  of  Time  for  Oblate  Spheroidal  Coordinates 


difference  is  required  throughout  the  PCM,  oblate  spheroidal  coordinates 
provide  a  better  choice  than  spheres.  The  effect  of  eccentricity  c/b  on  the 
time  required  for  complete  solidification  (TJ  =  0)  and  the  corresponding 
dimensionless  wall  temperature  are  qualitatively  similar  to  the  calculated 
results  for  the  elliptic  cylinder  case  shown  in  Figs.  6  and  7. 

For  phase-change  problems  in  an  outward  direction  where  7)  >  rjQ,  it 
can  be  shown  that  the  dimensionless  time  to  reach  a  given  position  rj  is 
obtained  by  placing  a  negative  sign  before  the  first  term,  inside  the  bracket, 
Eq.  (47).  The  wall  temperature  of  the  inner  boundary  T  can  be  calcu¬ 
lated  from  Eq.  (48)  by  placing  a  negative  sign  in  front  of  the  second  term 
inside  the  brackets.  The  resulting  expression  can  be  used  to  estimate  the 
rate  at  which  a  phase-change  point  would  move  into  the  half- space  above  the 
x-axis  from  a  circular  disc  of  radius  a,  located  on  the  x-axis,  which  is 
assumed  to  be  insulated  everywhere  else. 

E.  PROLATE  SPHEROIDAL  COORDINATES  (r],  Q,lJ/) 

The  prolate  spheroidal  coordinate  system  is  generated  by  rotating  an 
orthogonal  family  of  confocal  ellipses  and  hyperbolas  about  the  major  axis 
of  the  ellipses  as  shown  in  Fig.  2.  The  prolate  spheroids,  rj  =  constant  solids, 
vary  in  shape  from  near-spherical  to  thin  rods  of  finite  length.  The  rela- 

1  ^ 

tions  between  prolate  spheroidal  coordinates  and  Cartesian  coordinates  are 

x  =  a  sinh/;  sin#  cosl/' 
y  =  a  sinhTj  sin#  sin'l' 

z  =  a  cosh r]  cos#  (52) 

The  metric  coefficients  can  be  easily  derived  by  using  Eq.  (11) 

2  2  2 

=  gg  =  a  (sin  rj  +  sin  #) 

2  2  2 

g^  =  a  sin  T)  sin  #  (53) 
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The  thermal  resistance  between  two  prolate  spheroids  can  be  obtained  by 
using  Eq.  (22)  or  by  using  the  result  of  Ref.  15  directly. 


/n[tanh(77Q/2)  -  /n[tanh(7?/2)] 
47Tak 


(54) 


where  r)  <  T]q  for  the  inward  solidification  case.  The  following  derivation 
proceeds  in  a  similar  manner  to  the  oblate  spheroid  case,  namely,  writing 
the  heat  balance  equation  and  substituting  the  expressions  for  Aq  and  V 
which  are  functions  of  TJ  only.  The  dimensionless  form  of  the  heat  balance 
equation  is 


where  1^  and  are  two  integrals  that  can  be  integrated  in  closed  form  as 
shown  below: 


'o 

J *  ^sinh^rj  +  j  sinhr^dr;  = 


(cosh37)^  -  cosh3»7  -  cosh7)o  +  cosh rj) 


*0 

=  J*  /n[tanh(>7/2)]  sinh^77  +  y  sinh77)d77  = 


=  Y2  j^n[tanh(t)/2)]  (cos3/7q  -  cosh^) 

-  Y2  /n[tanh(T7/2)]  (cosh3?7  -  cosh rf) 

-  -J2  (cosh2r7Q  -  cosh2 T])  (57) 
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The  term  I  (77  )  is  a  definite  integral  required  to  define  the  surface  area  of  a 
prolate  spheroid.  Its  reciprocal  is  plotted  in  Fig.  3  as  a  function  of  77.  The 
Biot  number  is  defined  in  terms  of  the  semiminor  axis  c  which  is  related  to  a 

,  j'z  2 

by  a  =  v  b  -  c  . 

V ;  a  wall  temperature  can  be  obtained  in  a  similar  manner  to  the  oblate 
spheroiu  case;  the  result  is 

-1 


T  -  T 
w _ o 

T,  -  T 
fr 


1  +  Bi  sinh  TJ 


wf  - 1  • 


tanh(77  /2) 

/  _ o 

ln  tanh(J7/2) 


(58) 


Figure  9  shows  Eqs.  (55)  and  (58)  plotted  for  prolate  spheroids  of  c/b  ratios 
of  0.  9,  0.5,  and  0.  1.  The  frozen  fraction  F  is  related  to  the  77 -coordinate  by 


F  = 


V  -  V 
o 


l  0  sinh  TJcoshr] 

sinh2»7  coshT7 
o  o 


(59) 


Figure  9  illustrates  that  as  c/b  approaches  unity,  the  frozen  fraction  and 
dimensionless  wall  temperature  approach  the  limiting  case  of  a  sphere.  For 
small  values  of  c/b  the  limiting  results  are  close  to  the  c/b  =  0.  1  case  and 
cannot  be  accurately  estimated  by  the  infinite  cylinder  results.  For  cases 
where  the  temperature  difference  between  the  wall  and  the  freezing  point 
must  be  kept  at  a  low  value,  such  as  in  thermal  energy  storage  and  close 
temperature  control  applications,  the  phase-change  container  shapes  shown 
in  Fig.  9  are  not  desirable  because  they  all  result  in  the  maximum  possible 
temperature  difference  at  the  end  of  the  phase-change  process.  Elliptic 
cylinders,  oblate  spheroids,  and  other  slablike  geometries  result  n  lower 
temperature  differences  for  all  values  of  the  Biot  number. 

F.  BICYLINDRICAL  COORDINATES  (77,^,  z) 

The  bicylindrical  coordinate  system  is  generated  by  translating  the  two 
families  of  orthogonal  circles  in  the  xy-plane,  shown  in  Fig.  10,  parallel  to 
the  z-axis  pointing  into  the  paper.  The  circles  T)  =  const  are  drawn  about  two 
poles  x  =  ±a,  while  the  orthogonal  family  &  =  const  has  its  centers  on 
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+* 
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the  y-axis.  The  relation  between  bipolar  coordinates  and  Cartesian 
coordinates  is  ^ 


x  =  a  sinhTj/lposhTJ  -  cos  (A) 
y  =  a  sin'/'/fcoshT)  -  cos*/') 

z  =  z  (60) 


The  metric  coefficients  can  be  easily  derived  by  the  use  of  Eq.  (11) 

2 
a 

§  vj  _  ^  lb  ~  i  2  ’ 

'  ™  (cosh TJ  -  cosl/') 


g  =1 
sz 


(61) 


Using  the  thermal  resistance  between  two  circles  of  the  TJ  -  const  family 
( T)q  -  T))/2kln  and  proceeding  as  before  results  in 


k(T^r  -  TQ)t  I  sinhT^tanhTJ^  TJ° tanh  TJ^ 


pLw 


2Bi 


tanh^/7 


|csch^77  -  csch ^TJ 


rjc  sch  TJ  -  TJ  csch  TJ  J  +  ctnh TJ  -  ctnhty 


(62) 


where  w  is  the  location  of  the  center  of  the  17-circle  on  the  x-axis  and  r  is 
o  o 

the  radius  of  the  circle.  It  can  thus  be  shown  that  wq  =  acoth77o  and 
rQ  =  a/  |sinh77o|*  Equation  (62)  can  be  used  to  estimate  the  rate  of  phase 
change  between  eccentric  isothermal  cylinders  that  have  a  large  length-to- 
diameter  ratio.  As  TJ  approaches  zero,  Eq.  (62)  estimates  the  rate  of  phase 
change  around  a  long  pipe  of  radius  rQ  which  is  buried  at  a  depth  of  wq 
beneath  an  isothermal  level  surface  which  is  maintained  at  the  phase  change 
temperature. 
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G.  ESTIMATION  OF  ERRORS 

An  estimate  of  the  error  caused  by  neglecting  the  heat  capacity  of  the 
solidified  (or  melted)  portion  cannot  be  readily  determined  for  the  case  of 
the  elliptic  cylinder  and  spheroidal  coordinates.  However,  since  the  solu¬ 
tions  obtained  are  bracketed  by,  and  similar  to,  the  more  common  cases  of 
the  one -dimensional  slab  and  the  cylinder,  it  is  reasonable  to  assume  that  an 

analysis  of  the  error  involved  in  these  cases  will  indicate  the  error  involved 

2 

in  the  more  complex  geometries  treated  herein.  Goodman,  using  the  heat- 

balance  integral,  obtained  a  solution  for  the  slab  problem  shown  in  Figs,  lb 

and  Id  with  h  ~  oo  and  x  =0.  His  result  for  the  location  of  the  phase 
o  o  r 

change  front  is 

x  =  (yTTTlstTTT)1/2  y/zott  (63) 

which  approximates  the  exact  solution  given  in  Carslaw  and  Jaeger*  within 
2.  0%  up  to  Ste  =  0.  5.  The  comparable  solution  for  the  case  of  negligible 
heat  capacity  can  be  obtained  by  letting  hQ  approach  infinity,  taking  xq  =  0  in 
Eq.  (8),  and  carrying  out  the  simple  integration.  The  result  is 

x  =  >/  Ste  yiaT  (64) 

The  percent  error  in  calculating  the  time  required  to  reach  a  given  location 
of  the  phase -change  front  can  be  obtained  by  comparing  Eqs.  (63)  and  (64) 
(see  Fig.  11).  Ignoring  the  heat  capacity  results  in  change  of  phase  times 
that  are  always  shorter,  as  expected.  As  long  as  the  Stefan  number  is  below 
0.2,  Eq.  (64)  underpredicts  the  phase-change  time  by  less  than  10%.  For 
Stefan  numbers  below  0.  05  the  maximum  error  is  below  3%. 

Further  confirmation  of  the  validity  of  neglecting  the  sensible  heat 
effects  for  low  Stefan  numbers  was  carried  out  by  comparing  Megerlin's^ 
result  for  a  cylinder,  given  below,  with  Eqs.  (38)  and  (39). 
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ERROR  1%) 


Ste 


Fig.  11.  Percent  Error  in  SteFo  due  to  Neglecting  the  Sensible 
Heat  Effect  as  a  Function  of  Stefan  Number  with 
Biot  Number  as  a  Parameter 
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*  / - * - - - - 

SteFo^Q  =  J"  y(l  -  Bifnr  -  2  SteBi/nr  (2  -  Bi/nr  )  +  (1  -  Bi/nr  ) 


T  -  T 
w _ o_ 

T  -  T 
fr  o 


^  [r'rinr  +  ^(r'V)  (/nr)  ] 


dr 

(65) 

(66) 


where 


■  2SteBi 


r"  - 


2 

Bi^nr  )  -  2SteBi/nr  (2  -  Bi/nr  )  +  (1  -  B 


i/nr  ) 


These  equations,  although  approximate,  have  been  shown  to  yield  results  of 
high  accuracy  by  comparison  with  finite  difference  solutions, ^  The  percent 
error  as  a  function  of  Stefan  number  has  been  calculated  and  is  shown  in 
Fig.  11  with  Biot  number  as  a  parameter.  The  percent  error  is  defined  by 
100(SteFo(_.^Q  -  SteFo^._Q)/SteFo(_.^Q.  Lower  Biot  numbers  tend  to  decrease 
the  error  in  neglecting  the  sensible  heat  effects.  The  integration  of  Eq.  (65) 
was  carried  out  numerically  using  Simpson's  method  with  interval-halving 
and  a  convergence  criterion  of  0.  5  x  10  ^  between  iterations.  The  error  in 
comparing  the  wall  temperature  as  predicted  by  Eqs.  (66)  and  (39)  is  about 
the  same  as  above  for  Stefan  numbers  up  to  1.0.  For  larger  Stefan  numbers 
the  error  is  less  than  the  error  in  SteFo. 


H. 


FLUX  BOUNDARY  CONDITION 


When  the  heat  transfer  rate  is  specified  at  the  container  wall  the  heat 
balance  Eq.  (23)  simplifies  considerably  to 


Q  =  "Ldt(Vo  -  V| 


(67) 


Shamsundar,  N.,  and  Sparrow,  E.M.,  "Storage  of  Thermal  Energy  by 
Solid-Liquid  Phase-Change  - -Ten  perature  Drop  and  Heat  Flux,  11  Journal  of 
Heat  Transfer,  Trans.  ASME,  Series  C,  Vol.  96,  pp.  541-543. 
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for  the  inward  phase  change  problem  subject  to  the  same  assumptions  stated 
earlier.  The  wall  temperature  can  also  be  readily  calculated  from 
Q  =  (T^r  -  T^j/R  by  substituting  q^A^  for  Q  and  using  the  appropriate 
thermal  resistance  for  the  coordinate  system  used.  The  results  are  sum¬ 
marized  in  Table  1.  The  symbol  c  represents  the  half  thickness  xq  for  the 
case  of  a  slab,  the  outer  radius  r^  for  the  cylinder  and  sphere,  and  the 
semiminor  axis  for  the  remaining  coordinate  systems.  The  symbol  F  is 
related  to  T)  by  the  same  equations  as  for  the  convection  boundary  condition 
cases  described  earlier.  Figure  12  shows  the  calculated  results  for 
selected  geometric  shapes.  The  frozen  fraction  is  a  linear  function  of  the 
dimensionless  time  parameter  used  for  all  cases;  the  dimensionless  tempera 
ture  difference,  however,  is  nonlinear  and  tends  to  high  values  for  some 
geometries  at  the  end  of  the  phase  change  period.  For  applications  where 
low  temperature  differences  are  required,  elliptic  cylinders,  oblate 
spheroids,  and  slabs  are  more  appropriate  phase  change  container  shapes 
than  spheres,  cylinders  and  prolate  spheroids.  The  reason  is  that  the  ther¬ 
mal  resistance  between  the  container  wall  and  its  center  tends  to  infinity  for 
certain  geometries,  thus  resulting  in  high  temperature  differences. 

Note  that  the  equations  presented  in  Table  1  are  derived  from  thermal 
resistances  based  on  isothermal  boundary  conditions.  For  uniform  flux 
boundary  conditions  the  problem  becomes  generally  two-dimensional;  how¬ 
ever,  departures  from  one -dimensionality  are  expected  to  be  small  because 
of  the  symmetry  of  the  problems  considered  herein.  The  symbol  q  should  be 
interpreted  as  the  area-averaged  heat  transfer  rate  so  that  the  boundary 
condition  temperature  is  uniform.  This  is  believed  to  be  a  valid  approach 
since  in  most  real  situations  the  actual  boundary  condition  is  not  known 
exactly  and  usually  lies  somewhere  between  the  constant  temperature  and  the 
constant  flux  condition.  For  the  case  of  bipolar  coordinates  where  highly 
asymmetric  geometries  are  involved,  large  departures  from  one -dimensiona 
are  present  for  certain  uniform  flux  cases  as  reported  in  Ref.  17. 

^Thiyagarajan,  R.,  and  Yovanovich,  M.M.,  "Thermal  Resistance  of  a 
Buried  Cylinder  with  Constant  Flux  Boundary  Condition,  "  Journal  of 
Heat  Transfer,  Trans.  ASME,  Series  C,  Vol.  96,  pp.  249-250. 


Fig.  12.  Frozen  Fraction  and  Dimensionless  Temperature 
Difference  as  a  Function  of  Time  for 
Several  Geometries 


Table  1.  Flux  Boundary  Condition  Equations 


I.  APPROXIMATION  OF  TWO-DIMENSIONAL  PROBLEMS 

Sometimes  it  is  possible  to  approximate  the  thermal  behavior  of  two- 
and  three-dimensional  phase  change  problems  by  a  one-dimensional  analysis 
using  an  appropriate  coordinate  system.  Table  2  lists  a  comparison  of  the 
wall  temperature  ratio  and  dimensionless  time  parameter  for  different  frozen 
fractions  as  calculated  by  Eqs.  (38)  and  (39)  with  the  numerical  results 
obtained  in  Ref.  12  where  an  ad  hoc  computer  program  was  used  to  predict 
the  rates  of  freezing  and  temperature  distributions  of  a  PCM  contained  in  a 
convectively  cooled  square  container.  The  maximum  error  is  10%  for  this 
case  where  Biot  =  1.0.  The  maximum  error  calculated  for  the  Biot  =  10.0 
and  Biot  =  0.  1  cases  are  17  and  2.6%,  respectively.  Another  example  is  the 
two-dimensional  problem  of  a  long  cylinder  of  finite  length  that  can  be 
approximated  by  using  prolate  spheroidal  coordinates  where  c/b  approaches 
zero. 
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Table  2.  Comparison  of  Eqs.  (38)  and  (39)  with  the 
Numerical  Results  of  Ref.  12  for  Biot  =1.0 


Frozen 

Fraction 

(T  -  T  )/ (T ,  -T  ) 
w  o  '  fr  o 

SteFo 

Eq.  (3  9) 

Ref.  12 

Eq.  (38) 

Ref.  12 

0.  097 

0.9515 

0.95 

0. 04972 

0.05 

0.  189 

0.  9052 

0.  905 

0. 0993 

0.  10 

0.  274 

0.  8620 

0.  865 

0.  1474 

0.  15 

0.428 

0.  7817 

0.  815 

0. 2410 

0.25 

0.  564 

0.  7067 

0.  75 

0. 3325 

0.35 

0.  738 

0.5989 

0.64 

0.4658 

0.50 

0.  878 

0.4874 

0.  51 

0. 5940 

0.65 

0.  983 

0. 3292 

0.  36 

0.7199 

0.  80 

III.  DISCUSSION  AND  CONCLUSIONS 


A  general  expression  for  predicting  the  rates  of  melting  or  freezing 
and  the  resulting  temperature  histories  was  presented,  based  on  the  conduc¬ 
tion  shape  factor  equation  derived  in  Refs.  13  or  15  and  the  assumptions 

7 

adopted  by  London  and  Seban  ,  concerning  the  phase  change  process. 
Closed-form  solutions  were  obtained  for  elliptic  cylinder,  oblate  spheroidal, 
prolate  spheroidal  coordinates,  and  bicylindrical  coordinates.  Arithmetic 
results  were  presented  in  graphical  form  for  the  inward  phase  change  prob¬ 
lem  in  some  of  these  coordinate  systems.  These  results  are  of  special 
interest  to  the  design  of  thermal  storage  equipment  where  container  wall 
temperatures  and  heat  flow  rates  for  various  container  shapes  are  important 
as  is  the  location  of  the  solid-liquid  interface.  Elliptic  cylinders  and  oblate 
spheroids  of  medium  to  high  eccentricities  that  resemble  a  slab-like  geo¬ 
metry  were  shown  to  give  lower  temperature  differences  and  more  uniform 
heat  flow  rates  than  cylinders,  spheres,  and  prolate  spheroids.  Low  Biot 
numbers  resulted  in  more  uniform  heat  extraction  rates  and  low  temperature 
differences;  however,  longer  times  are  required  for  complete  phase  change. 
The  error  involved  in  neglecting  the  sensible  heat  effects  was  computed  in 
order  to  assess  the  limits  of  applicability  of  the  results  presented.  The 
maximum  expected  error  was  determined  to  be  less  than  10%  for  Stefan 
numbers  below  0.2  and  3%  for  Stefan  numbers  less  than  or  equal  to  0.05. 

The  large  number  of  other  applications  where  the  solution  technique 
can  be  used  is  of  far  greater  significance  than  the  few  examples  calculated 
herein.  For  example,  freezing  around  a  disc-shaped  cryoprobe  can  be  esti¬ 
mated  by  using  the  expressions  presented  in  oblate  spheroidal  coordinates  in 
an  outward  direction.  Freezing  or  melting  around  a  finite  plate  immersed 
in  a  semi-infinite  medium  can  be  estimated  by  the  elliptic  cylinder  results 
in  an  outward  direction.  Other  cases  such  as  melting  or  freezing  around 
buried  cylinders  or  spheres  can  be  treated  by  the  use  of  bicylinder  and  bi- 
spherical  coordinate  systems,  respectively. 
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NOMENCLATURE 


A 

a 

Bi 

b 

C 

c 

F 

Fo 

g 

h 

k 

L 

Q 

q 

R 

r 

V 

r 

r* 

Ste 

SteFo 

T 

t 


surface  area  of  container 
focal  length  of  ellipse 
Biot  number,  hc/k 
semimajor  axis 
specific  heat 

semiminor  axis  and  specific  heat 
frozen  fraction 

2 

Fourier  number,  (k/pCx^jt 
metric  coefficient 

convective  heat  transfer  coefficient 

thermal  conductivity 

latent  heat  of  fusion 

total  heat  transfer  rate 

heat  transfer  rate  per  unit  area 

thermal  resistance 

radius  of  circular  cylinder  or  sphere 
dimensionless  radius  for'cylinder  or  sphere 
time  rate  of  change  of  solid- liquid  interface  radius 
Stefan  number,  C(Tf  -  T  )/L 

Product  of  Stefan  and  Fourier  numbers,  k(T^r  -  TQ)t/pLr 

temperature 

time 
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NOMENCLATURE  (Continued) 


ul’  U2*  U3 
x.y,  z 

a 

r\,e,z 

*1,0,* 

p 

i 

Subscripts 

c 

fr 

o 


orthogonal  curvilinear  coordinates 

Cartesian  coordinates 

thermal  diffusivity,  k /PC 

elliptic  cylinder  coordinates 

oblate  or  prolate  spheroidal  coordinates 

density 

length 


complete  freezing  or  melting 

corresponding  to  the  coordinate  i,  j,  or  k 

freezing  (or  melting)  temperature 

outer  wall  for  inward  problems 
inner  wall  for  outward  problems 
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